Pitchfork bifurcations in blood-cell shaped dipolar Bose-Einstein condensates 
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We demonstrate that the method of coupled Gaussian wave packets is a full-fledged alternative to 
direct numerical solutions of the Gross-Pitaevskii equation of condensates with electromagnetically 
induced attractive 1/r interaction, or with dipole-dipole interaction. Moreover, Gaussian wave 
packets are superior in that they are capable of producing both stable and unstable stationary 
solutions, and thus of giving access to yet unexplored regions of the space of solutions of the Gross- 
Pitaevskii equation. We apply the method to clarify the theoretical nature of the collapse mechanism 
of blood-cell shaped dipolar condensates: On the route to collapse the condensate passes through 
a pitchfork bifurcation, where the ground state itself turns unstable, before it finally vanishes in a 
tangent bifurcation. 
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Bose-Einstein condensates (BECs) with dipole-dipole 
interaction have become an active and exciting field of re- 
search because they offer the possibility of tuning the rel- 
ative strengths of the short-range isotropic contact inter- 
action and the anisotropic long-range dipole interaction 
by manipulating the s-wave scattering length via Fesh- 
bach resonances, and thus of studying a wealth of new 
phenomena that occur as one crosses the whole range 
from dominance of the contact interaction to that of 
the dipole interaction. The experimental realization of 
a BEC of chromium atoms [IH3], which possess a strong 
magnetic dipole moment, has given additional impetus 
to the field (for a comprehensive list of references see 
the recent review by Lahaye et al. |4j). In the dilute 
limit, the theoretical description of these condensates can 
be done in the framework of the Gross-Pitaevskii equa- 
tion (GPE). This nonlinear Schrodinger equation has 
been solved in the literature so far by simple variational 
ansatzes, where the mean-field energy is minimized, e.g., 
with the conjugate gradient method or by imaginary time 
evolution. In this Letter we will show that the method 
of coupled Gaussian wave packets is an adequate alterna- 
tive to solving the GPE of BECs with long-range interac- 
tions. Moreover, we will show that the method is superior 
in that it also yields unstable stationary solutions, and 
thus opens access to regions of the space of solutions of 
the GPE unexplored heretofore. As an application of the 
method we will analyze in detail the theoretical nature 
of the collapse mechanism of dipolar BECs. 

The GPE for ultracold gases with long-range interac- 
tions, described by the interatomic potential Vi r (r), has 
the form 

Z rJ^ (r ' = [ ~ A + + T ' y2 + 

+8TrNa\ijj(r,t)\ 2 + V h (r)]ijj(r,t) , (1) 

where for dipolar interaction we have 
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with i?' the angle between r — r' and the axis of an exter- 
nal magnetic field. For completeness we will also consider 
the case of an isotropic "gravity-like" attractive 1 jr long- 
range interaction, 

W .- 2 WdV^|. (3) 
J \r — r'\ 

According to O'Dell et al. [5] this interaction could be 
electromagnetically induced by exposing the condensate 
atoms to an appropriately arranged set of triads of laser 
beams. The appealing feature of such "monopolar" 
condensates is that they can be self-trapping, i.e. ex- 
ist without an external trapping potential. The equa- 
tions above have been brought into dimensionless form 
by introducing natural units, which for monopolar in- 
teraction (V m0 no = —u/r) are [5H5) the "Bohr radius" 
o-u = h 2 /{mu) for lengths, the "Rydberg energy" E u — 
h 2 /(2ma 2 l ) for energies and u u = E u /h for frequen- 
cies. Natural units for dipolar atoms with magnetic mo- 
ment n are [9] the dipole length ad = /xo/xm/(27rfi, 2 ), the 
dipole energy — h 2 /(2ma 2 i ) and the dipole frequency 
= Ed/h. The quantities Jx,y,z in denote the trap- 
ping frequencies in the three spatial directions measured 
in the respective frequency units, N is the number of 
bosons, and a the scattering length in units of a u and a<j, 
respectively. 

The most obvious way of solving the Gross-Pitaevskii 
equation is its direct numerical integration on multi- 
dimensional grids using, e.g., fast- Fourier techniques. 
The stationary ground state can be obtained by imagi- 
nary time evolution. These calculations, however, may 
turn out laborious, and physical insight can often be 
gained using approximate, in particular, variational solu- 
tions. A common approach employed for determining the 
dynamics and stability of condensates both with contact 
interaction only [101 lll| and with additional long-range 
interaction jlSj is to assume a simple Gaussian form of 
the wave function, with time-dependent width parame- 
ters and phase, and to investigate the dynamics of these 
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quantities. For dipolar condensates improvements on the 
simple Gaussian form were made by multiplying it by 
second-order Hermite polynomials |12j . 

As an alternative to numerical quantum simulations 
on multidimensional grids we will extend the variational 
calculations in such a way that numerically converged 
results are obtained with significantly reduced computa- 
tional effort compared to the exact quantum simulations 
but with similar accuracy. The method is that of cou- 
pled Gaussian wave packets. It was originally proposed 
by Heller [131 114] to describe quantum dynamics in the 
semiclassical regime, and was successfully applied to the 
dynamics of molecules and atoms in external fields |15j . 
The idea is to choose trial wave functions which are super- 
positions of N different Gaussians centered at the origin 

N N 
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where both the width parameters a k and the scalars 7 fc 
are complex quantities, with the latter determining the 
weight and the phase of the individual Gaussian. In- 
serting the ansatz Q into the time-dependent Gross- 
Pitaevskii equation and applying the time-dependent 
variational principle where — Hip(t)\\ 2 is mini- 

mized by varying <E>, and afterwards $ is set equal to 
$ = V», yields a set of ordinary differential equations for 
the width parameters a k and the scalars (cf. [T5] ) 

4 = -4(^) 2 - -V 2 k ; (3 = x, y, z; (5a) 

V '=2t(o* + a* + o*)-t;*. (5b) 

The quantities (vq , V 2 fe ) with k = 1, . . . , N constitute the 
solution vector to the set of linear equations 
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with I — 1, . . . , N; m + n = 0, 2; and X\ = x, x 2 = y, 
x 3 = z. Here, V(x) = V c + Vh- + V t denotes the sum of the 
contact, the long-range and of external trap potentials. 
The important and appealing point of this procedure is 
that all necessary integrals with the trial wave functions 
g l , g k from Q can be calculated analytically. 
Stationary variational solutions to the extended Gross- 
Pitaevskii equation ([I]) are found by searching for the 
fixed points of jH]), i.e. solving a k — 0:/f k — for each 
k = 1 , . . . , N via a 4 N dimensional highly nonlinear root 
search. The resulting stationary width and weight /phase 
parameters can then be used to calculate the mean field 
energy E m{ = (*| - A + V t + §(V C + Vt)|*> and the 
chemical potential \i = (*| - A + V t + V c + V^). 



It is important to note that in contrast to numerical 
calculations with imaginary time evolution, which only 
work for stable solutions, this procedure will produce 
both stable and unstable solutions, and thus uncover 
yet unexplored parts of the space of solutions of the 
Gross-Pitaevskii equation. 

To analyze the stability of the stationary solutions the 
dynamical equations ^ are split into real (R) and imag- 
inary (/) parts and linearized around the fixed points. 
The eigenvalues of the Jacobian matrix J at the fixed 
point 

0(a^V/,y«y.') ' (7j 

with a, (3 = x, y, z; k,l — 1, N, determine the stability 
properties of the solution. If all eigenvalues Xj of the sys- 
tem are purely imaginary, the motion is confined to the 
vicinity of the fixed point and quasi-periodic. If one real 
part or several real parts of the eigenvalues are non-zero, 
small variations from the fixed point lead to exponential 
growth of the perturbation. 
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FIG. 1: Chemical potential /i for self-trapped condensates 
with attractive 1/r interaction as a function of the scaled 
scattering length N 2 a obtained by using up to 5 Gaussian 
wave packets in comparison with the result of the exact nu- 
merical solution of the stationary Gross-Pitaevskii equation. 
Note that all forms yield a tangential bifurcation diagram, 
with a stable (upper) and an unstable (lower) branch. The 
inclusion of three Gaussians already well reproduces the exact 
numerical result. 

As a first application we demonstrate the efficiency of 
the coupled Gaussian wave packet method for conden- 
sates with attractive 1/r long-range interaction, for the 
case of self-trapping {^ x ,y,z = 0). Figure [l] shows the 
results for the chemical potential as a function of the 
scaled scattering length N 2 a for superpositions of 1 to 
5 Gaussians in comparison with the results of the exact 
numerical solution. It is evident that all forms reproduce 
the bifurcation behavior discussed in [5HH] : at a critical 
point two solutions of the Gross-Pitaevskii equation are 
born in a tangent bifurcation, one stable (upper branch) 
and one unstable (lower branch). The numerically ac- 
curate bifurcation point lies at a CT w —1.025147. It can 
also be seen that, while the variational calculation with 
one Gaussian, with a£| =1 = —1.178, still lies far off the 
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correct result, the inclusion of only one more Gaussian 
brings the chemical potential curve already close to the 
numerical result, and practically no improvement is vis- 
ible in Fig. [T] when 3 or more Gaussians are included. 
Using 5 coupled Gaussians the exact bifurcation point is 
reproduced with an accuracy of 10~ 6 . Similar results are 
obtained in the presence of a trapping potential. 

We now turn to dipolar condensates. Previous studies 
[T2l ITS] have shown that in certain regions of the param- 
eter space dipolar condensates assume a non-Gaussian 
biconcave "blood-cell-like" shape. To demonstrate the 
power of the coupled Gaussian wave packet method, we 
choose a set of such parameters. We consider an axisym- 
metric trap with (particle number scaled) trap frequen- 
cies N 2 j z = 25200 along the polarization direction of 
the dipoles and N 2 j p = 3600 in the plane perpendicular 
to it (corresponding to an aspect ratio of A = "fz/jp = 
7). For this set of parameters we show in Fig. [2] (a) the 
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FIG. 2: (a) Convergence of the mean-field energy with in- 
creasing number of coupled Gaussian wave packets (squares) 
and comparison with the value obtained by a lattice calcula- 
tion with grid size 128 x 512 (dashed line), which lies ener- 
getically higher than the exact converged variational solution 
(solid line), (b) Comparison of the variational wave func- 
tion for 6 coupled Gaussians (solid curves) with values of the 
numerical one (triangles) at different z coordinates. Both so- 
lutions show a biconcave shaped condensate. The figures are 
for (particle number scaled) trap frequencies N 2r y z = 25200 
and N 2 j p = 3600, and scattering length a = 0. 



convergence behavior of the mean field energy. We com- 
pare the variational solution as the number of Gaussian 
wave packets is increased from 2 to 6 with the mean field 
energy value of a numerical lattice calculation (imagi- 
nary time evolution combined with FFT) with a grid 
size of 128 x 512, at scattering length a = as an 
example. The mean field energy for one Gaussian is 
E m { — 60361 Ed and lies far outside the vertical energy 
scale. Evidently the numerical value is more than excel- 
lently reproduced by 5 and 6 coupled Gaussians. The 
behavior for other scattering lengths is similar. Also the 
wave function nicely converges, and moreover, as can be 



seen in Fig. [2] (b), reproduces the biconcave shape of 
the condensate as does the numerical solution. Thus the 
method of coupled Gaussians is a viable and full-fledged 
alternative to direct numerical solutions of the Gross- 
Pitaevskii equation for dipolar condensates. 
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FIG. 3: (a) Mean field energy of a dipolar condensate for 
(particle number scaled) trap frequencies N 2r y z — 25200 and 
N 2r y p = 3600 as a function of the scattering length. In the 
variational calculation with one Gaussian a stable ground 
state (g N ~ X ) and an unstable excited state (e N=1 ) emerge in 
a tangent bifurcation. Using coupled Gaussians two unstable 
states emerge (labeled u coup ed ), Q f which the lower one turns 
into a stable ground state (g cou P led ) i n a pitchfork bifurca- 
tion, (b) Stability eigenvalues A of the pitchfork bifurcation 
point for calculations with 6 coupled Gaussians, scattering 
length in rectangle marked in (a). Real and imaginary parts 
of two selected eigenvalues of the Jacobian (J7| as a function 
of the scattering length. For a < aj? r = —0.00359 the solu- 
tion is unstable with one pair of real eigenvalues. At a P r the 
real eigenvalues vanish in a pitchfork bifurcation and a stable 
ground state forms with purely imaginary eigenvalues. Only 
those eigenvalues involved in the stability change are shown. 

Figure [3] (a) shows, for the same set of trap frequen- 
cies, the results for the mean field energy of the conden- 
sate as a function of the scattering length a (in units ad) 
for a wave function with one Gaussian, and for 5 cou- 
pled Gaussian wave packets. Results obtained using 6 
Gaussians would be indistinguishable in the figure from 
those obtained using 5 Gaussians, and the results for 2- 
4 Gaussians are not shown for the sake of clarity of the 
figure. 

Similar to the above findings for monopolar conden- 
sates, and as is known from previous variational calcula- 
tions [§] for dipolar condensates, for N = 1 two branches 
of solution are born in a tangent bifurcation. The en- 
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ergetically higher branch has purely real stability eigen- 
values ±A fl , corresponding to an unstable excited state 
e N=1 , the lower branch possesses purely imaginary eigen- 
values ±A J and corresponds to the stable ground state 
g N=1 . At the bifurcation point the branches of the sta- 
bility eigenvalues merge and vanish. 

The situation is different if the condensate wave func- 
tion is described by more than one Gaussian. As the 
scattering length is decreased from positive values to- 
wards the tangent bifurcation, the branch corresponding 
to the ground state g cou P led turns into an unstable state 
u cou P ied at a scatter i n g length of a? = -0.00359. This 
is evident from the stability analysis shown in Fig. [3] (b) 
where the stability eigenvalues for the ground state, cal- 
culated using 6 Gaussians, are plotted in a small interval 
of the scattering length around a% T . Above a P r the eigen- 
values are purely imaginary, below they are purely real. 
[Note that in a Bogoliubov analysis this instability should 
appear as a dynamical instability] The ground state re- 
mains unstable down to the tangent bifurcation point at 
a* r = —0.01224, where it joins the branch of the unstable 
excited state. 

The quality of the calculation using 5 Gaussian wave 
packets is also demonstrated in Fig. [3] (a) where the re- 
sults of a numerically grid calculation by imaginary time 
evolution are shown by crosses. Evidently the numerical 
results and the results obtained using 5 coupled Gaus- 
sians excellently agree. The imaginary time calculation, 
however, can only trace the stable branch of the solution 
and fails for the unstable branch. Thus it is demonstrated 
that the Gaussian wave packet method is not only numer- 
ically accurate but also capable of giving access to regions 
of the space of solutions of the Gross-Pitaevskii equation 
with dipolar interaction that are difficult to investigate 
by conventional numerical full quantum calculations. 

The phenomenon of one smooth branch of solutions 
becoming unstable as a function of a control parame- 
ter is reminiscent of a pitchfork bifurcation. The two 
stable solutions on the fork arms which should also be 
born in a pitchfork bifurcation, and exist in a tiny neigh- 
borhood (a P r — e) < a < a p r , are numerically hard to 
trace and therefore not plotted in the figure. Their ex- 
istence, and the pitchfork type of the bifurcation, how- 
ever, can be proven by looking at the "phase portrait" 
plotted in Fig. [4] at a value of the scattering length 
a = —0.036 slightly below a p r . Figure [4] shows con- 
tours of equal deviation of the mean field energy from 
that of the ground state in the plane spanned by the two 
eigenvectors whose eigenvalues are involved in the sta- 
bility change in Fig [3] (b). The coordinate axes 61,62 
correspond to small variations of the Gaussian parame- 
ters in the eigenvector directions around the hyperbolic 
fixed point solution located at the origin. The portrait 
clearly reveals the existence of two nearby elliptic fixed 
points corresponding to two additional stable solutions. 
Therefore, in a small interval e of the scattering length 



below a p ., (a P r — e) < a < a P r , there exist two additional 
branches, besides the unstable solution. This proves that 
the bifurcation is of pitchfork type. Note that the clas- 
sification of the condensate as unstable for a < a P r nev- 
ertheless remains true in physical terms due to the nu- 
merically small value of e. We also note that for a > a P r 
the phase portrait possesses only one elliptic fixed point 
cooresponding to the stable stationary ground state. 




FIG. 4: Contour plot of the mean field energy with the eigen- 
vectors corresponding to the eigenvalues of Fig.[3](b) lineariz- 
ing the vicinity of the fixed point (Si , 82 in arbitrary units) . 
The figure shows a — —0.0036 close below the pitchfork bi- 
furcation point, showing three fixed points: Two stable and 
one hyperbolic. 

Is there a chance of observing dipolar BECs on the 
stable fork arms? The answer probably is no, in the 
same way as it is in the case of the question of observing 
the transition to structured ground states, possibly as- 
sociated with a roton instability, shortly before collapse. 
The reason is the difficulty of adjusting trap frequencies 
and the scattering length to the necessary precision in a 
real experiment. Nevertheless theoretical investigations 
of this type close to the threshold of instability of dipo- 
lar condensates are valuable in their own right since they 
help to understand the nature of the collapse, and thus 
of "what's really going on" . 
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